Dimension reduction and representation learning

Adapted from Childhood Cancer Data Lab, Alex’s Lemonade Stand Foundation training materials.

Background

Before we get into it, let’s briefly cover some concepts.

Dimension reduction (also called dimensionality reduction) is any technique that reduces the number of features or dimension in a dataset. A dimension in our high-dimensional space is the measurement of an individual gene in our RNA-seq data in our particular case. Dimension reduction can come about through feature elimination or feature selection, where we remove features or variables if they are unlikely to give us meaningful information (or are redundant with other variables). Or we can use feature extraction, where we form new features from the input data. (See Freytag (2019) and Raj (2019).)

In general, the idea behind dimension reduction, as we’ll cover it in this notebook, is that we want to transform our input data into a low-dimensional representation that retains some of the meaningful structure (feature extraction).

If we want to get the sense of the overall structure in our dataset, sometimes we want the following:

Representation learning (also called feature engineering or feature learning) is any technique that automatically detects representations from input data. We’ll cover unsupervised methods like matrix factorization here, but supervised methods like supervised neural networks are examples of representation learning, too.

When we’re working with gene expression data, often what we want is the following:

We may want to use these representations in a semi-supervised manner, where we learn representations from a collection of data without labels and then use them to encode a low-dimensional representation of data with labels that we then use for prediction.

Set up

# Set seed for reproducibility
set.seed(2020)

Libraries

library(tidyverse)

Registered S3 methods overwritten by 'dbplyr':
  method         from
  print.tbl_lazy     
  print.tbl_sql      
── Attaching packages ─────────────────────────────────────────────────────────────────────────────── tidyverse 1.3.0 ──
✓ ggplot2 3.3.2     ✓ purrr   0.3.4
✓ tibble  3.0.3     ✓ dplyr   1.0.0
✓ tidyr   1.1.0     ✓ stringr 1.4.0
✓ readr   1.3.1     ✓ forcats 0.5.0
── Conflicts ────────────────────────────────────────────────────────────────────────────────── tidyverse_conflicts() ──
x dplyr::filter() masks stats::filter()
x dplyr::lag()    masks stats::lag()

# UMAP implementation
library(uwot)

Loading required package: Matrix

Attaching package: ‘Matrix’

The following objects are masked from ‘package:tidyr’:

    expand, pack, unpack

# t-SNE implementation
library(Rtsne)

Files

clinical_file <- file.path("data", "pbta-medullo-histologies.tsv.gz")

# Processed RNA-seq data
rnaseq_file <- file.path("data", "pbta-medullo-vst-collapsed.tsv.gz")

# Pathway-Level Information ExtractoR model we've prepared ahead of time
plier_file <- file.path("models", "pbta-medullo-plier.RDS")

Read in data

Read in our transformed RNA-seq data.

rnaseq_df <- read_tsv(rnaseq_file)

Parsed with column specification:
cols(
  .default = col_double(),
  gene_symbol = col_character()
)
See spec(...) for full column specifications.

For the applications in this notebook, we’ll transpose RNA-seq data.

rnaseq_mat <- rnaseq_df %>%
  column_to_rownames("gene_symbol") %>%
  as.matrix() %>%
  t()

Because we are interested in the overall structure of our data, we’re going to subset the matrix to features with high variance moving forward. Eliminating features that don’t vary within themselves and therefore are less likely to contribute to the overall structure in our data that we are most interested in. Fewer features is appealing from a computational resource standpoint!

# Genes are now columns
gene_variance <- matrixStats::colVars(rnaseq_mat)
# We'll want to retain the top 10%
variance_threshold <- quantile(gene_variance, 0.95)
# Filter to top 10%
high_var_mat <- rnaseq_mat[, gene_variance > variance_threshold]

Read in our clinical metadata file.

clinical_df <- read_tsv(clinical_file)

Parsed with column specification:
cols(
  Kids_First_Biospecimen_ID = col_character(),
  short_histology = col_character(),
  broad_histology = col_character(),
  tumor_descriptor = col_character(),
  molecular_subtype = col_character()
)

For the most part, we’ll only use the molecular subtype information from the clinical meatadata.

subtype_df <- clinical_df %>%
  select(Kids_First_Biospecimen_ID, molecular_subtype)

Dimension reduction

Adapted from CCDL scRNA-seq dimension reduction and Dimensionality Reduction concepts

Principal Components Analysis (PCA)

Principal components analysis (PCA) is a dimensionality reduction technique that captures the main sources of variation in our data in the first two principal components (PC1 and PC2). PCA applies linear algebra to find the axis of greatest variation in the data matrix, then the axis perpendicular to that with the most remaining variation, then the one perpendicular to those with the most remaining, etc.

# PCA
pca_results <- prcomp(high_var_mat, scale = TRUE)

# Join the first two PCs with the subtype labels for plotting purposes
pca_df <- data.frame(pca_results$x[, 1:2]) %>%
  rownames_to_column("Kids_First_Biospecimen_ID") %>%
  inner_join(subtype_df, by = "Kids_First_Biospecimen_ID")

# Create a plot of the first two PCs
ggplot(pca_df,
       aes(x = PC1, 
           y = PC2,
           color = molecular_subtype)) +
  geom_point(alpha = 0.75) +
  colorblindr::scale_color_OkabeIto() +
  theme_bw()

PCA can sometimes produce results where the points occupy the whole plotting space, which would likely be more pronounced if we had even more samples.

Uniform Manifold Approximation and Projection (UMAP)

UMAP (Uniform Manifold Approximation and Projection) is a machine learning technique designed to provide more detail in highly dimensional data than a typical PCA. While PCA assumes that the variation we care about has a particular distribution (normal, broadly speaking), UMAP allows more complicated distributions that it learns from the data. We won’t go into the underlying mathematics, but you can look at the paper by McInnes, Healy, & Melville (2018). The main advantage of this change in underlying assumptions is that UMAP can do a better job separating clusters, especially when some of those clusters may be more similar to each other than others.

Another dimensionality reduction technique that you may have heard of is t-SNE (t-distributed Stochastic Neighbor Embedding), which has similar properties to UMAP, and often produces similar results. There is some ongoing debate about which of these two techniques is superior for high-dimensional biological data, and whether the differences are due to the underlying algorithm or to implementation and parameter initialization defaults. Regardless of why, in our experience, UMAP seems to produce slightly better results in this context and run a bit faster, but the differences can be subtle.

Let’s use UMAP with the default parameters and see what the plot looks like!

# Run UMAP with default parameters, but perform PCA first and use the specified
# number of PCs. This is recommended for > 100 columns. 50 PCs are often 
# recommended for t-SNE applications, which we'll see in a moment.
umap_results <- uwot::umap(high_var_mat, pca = 50)

# Add the subtype labels to the coordinates returned by uwot::umap
umap_df <- data.frame(umap_results) %>%
  rename(UMAP1 = X1, UMAP2 = X2) %>%
  bind_cols(subtype_df)

# Make a plot that's very similar to the one we made for PC1 and PC2
ggplot(umap_df,
       aes(x = UMAP1, 
           y = UMAP2,
           color = molecular_subtype)) +
  geom_point(alpha = 0.75) +
  colorblindr::scale_color_OkabeIto() +
  theme_bw()

Now that we have an idea of what a UMAP plot with the default parameters looks like, let’s try experimenting with the n_neighbors parameter (number of nearest neighbors). This parameter constrains the size of the local neighborhood size when learning the manifold structure. Lower values emphasize the local structure of the data (ref).

In the interest of experimenting with the n_neighbors parameter, we’ve written a wrapper function that will allow us to

# For a given n_neighbor value, run UMAP and make a plot using the coordinates
# it returns with samples colored by their subtype label.
umap_plot_wrapper <- function(n_neighbors_param) {
  
  # Perform UMAP using the specified n_neighbors value
  umap_results <- uwot::umap(high_var_mat, 
                             n_neighbors = n_neighbors_param,
                             pca = 50)

  # Add the subtype labels
  umap_df <- data.frame(umap_results) %>%
    rename(UMAP1 = X1, UMAP2 = X2) %>%
    bind_cols(subtype_df)
  
  # Make a plot that's very similar to the one we made for PC1 and PC2
  ggplot(umap_df,
         aes(x = UMAP1, 
             y = UMAP2,
             color = molecular_subtype)) +
    geom_point(alpha = 0.75) +
    colorblindr::scale_color_OkabeIto() +
    theme_bw()
}

Let’s start by running UMAP using the same n_neighbors value but a different seed.

# Set a different seed but use the same n_neighbors
set.seed(1234)
umap_plot_wrapper(n_neighbors_param = 15)

For visualization, this difference doesn’t change our interpretation, but it’s a good idea to run this multiple times with different seeds to make sure you’re not getting anomalous results.

Okay, let’s try a variety of n_neighbor values and see what happens.

umap_plot_wrapper(n_neighbors_param = 3)

umap_plot_wrapper(n_neighbors_param = 30)

umap_plot_wrapper(n_neighbors_param = 50)

t-Distributed Stochastic Neighbor Embedding (t-SNE)

Like UMAP, t-SNE is nonlinear dimensionality technique that can produce well-separated clusters.

# This implementation performs PCA first.
# We are using the default perplexity value here, but like we saw with
# n_neighbors this parameter is highly influential on your results!
tsne_results <- Rtsne::Rtsne(high_var_mat, perplexity = 30)

# Add the subtype labels to the coordinates returned by Rtsne::Rtsne()
tsne_df <- data.frame(tsne_results$Y) %>%
  rename(tSNE1 = X1, tSNE2 = X2) %>%
  bind_cols(subtype_df)

# Make a plot using the coordinates from t-SNE
ggplot(tsne_df,
       aes(x = tSNE1, 
           y = tSNE2,
           color = molecular_subtype)) +
  geom_point(alpha = 0.75) +
  colorblindr::scale_color_OkabeIto() +
  theme_bw()

Some takeaways!

• Newer methods that learn nonlinear patterns can help us find well-separated clusters, which is handy for visualization.
• The results strongly depend on parameters that balance local structure vs. global structure, so these values should be chosen based on the question on hand and it may be useful to test a range of values.
• The interpretation of these methods may sometimes be counter-intuitive (e.g., distance between clusters).

A domain-specific dimension reduction technique: Pathway-Level Information ExtractoR (PLIER)

Adapted from CCDL OpenPBTA PLIER and Plotting latent variables notebooks.

All of the techniques we’ve seen so far can (can, not should!) be applied to any generic set of data. Now we will cover a method called Pathway-Level Information Extractor (PLIER) (Mao et al. (2019)).

Let’s revisit the example from above and put it in the context of PLIER:

We like PLIER for a few reasons:

• It is a matrix factorization approach. That means we can get a low-dimensional representation of our data. Specifically, PLIER learns correlated patterns of expression in our data or latent variables (LVs). Here, a latent variable is an “eigengene-like” combination of genes’ expression. (It’s called latent because it’s not directly measured, but instead inferred from the individual gene measurements.)
• It includes penalties such that some of the LVs will align with gene sets that we give it, so it’s excellent for biological discovery.
• The authors demonstrated that it performs favorably with regard to estimating proportion of immune cells in a sample as compared to other methods.
• Because not all LVs align with the gene sets we input, some of them capture unwanted technical variation. In our experience with the method, it does this quite well.

Here’s an overview of the PLIER method from Mao et al. (2019) (Figure 1).

Fig. 1 | PLIER overview. PLIER is a matrix factorization approach that decomposes gene expression data into a product of a small number of LVs and their corresponding gene associations or loadings, while constraining the loadings to align with the most relevant automatically selected subset of prior knowledge. a, Given two inputs, the gene expression matrix Y and the prior knowledge (represented as binary gene set membership in matrix C), the method returns the LVs (B), their loadings (Z), and an additional sparse matrix (U) that specifies which (if any) prior-information gene sets and pathways are used for each LV. The light gray area of U indicates the large number of zero elements of the matrix. We apply our method to a whole-blood human gene expression dataset. b, The positive entries of the resulting U matrix are visualized as a heat map, facilitating the identification of the correspondence between specific LVs and prior biological knowledge. As the absolute scale of the U matrix is arbitrary, each column is normalized to a maximum of 1. c, We validate the LVs mapped to specific leukocyte cell types by comparing PLIER estimated relative cell-type proportions with direct measurements by mass cytometry. Dashed lines represent 0.05, 0.01, and 0.001 significance levels for Spearman rank correlation (one-tailed test). NK cell, natural killer cell.

In the interest of time, we’ve prepared the model in advance that we will read in now. You can see the steps we took in setup/02-plier-prep.Rmd.

plier_results <- read_rds(plier_file)

PLIER gives us information about the association between a latent variable it learns and an input gene set in the summary data frame. In the interest of time, we’re going to focus in on a single latent variable to illustrate how you could use these values.

plier_results$summary %>% 
  filter(`LV index` == 56,
         FDR < 0.05) %>%
  arrange(FDR)

The B matrix contains the estimates or expression values for all of the latent variables learned by the model.

dim(plier_results$B)

[1]  58 121

Let’s get this into long format and add the subtype labels.

lv_df <- data.frame(plier_results$B) %>%
  # Data frame where the latent variable names are in a column
  rownames_to_column(var = "LV") %>%
  # Wide to long
  pivot_longer(cols = starts_with("BS_"),
               names_to = "Kids_First_Biospecimen_ID",
               values_to = "LV_estimate") %>%
  # Add the subtype labels
  inner_join(subtype_df, by = "Kids_First_Biospecimen_ID")

And plot the estimates for LV56 – we can use this to get an idea about the differences in immune-related signals between the subtypes.

lv_df %>%
  # PLIER names the latent variables based on the gene sets that they're
  # associated with, but as we saw above this is not the only gene set that's
  # associated with this LV!
  filter(LV == "56,SVM Macrophages M2") %>%
  # Make a boxplot where samples are grouped by their molecular subtype and 
  # that shows the individual samples as points
  ggplot(aes(x = molecular_subtype, 
             y = LV_estimate, 
             group = molecular_subtype, 
             colour = molecular_subtype)) +
  geom_boxplot(outlier.shape = NA) +
  geom_jitter(width = 0.2, alpha = 0.5) +
  # Use the same theme and color scheme we have been using throughout
  theme_bw() +
  colorblindr::scale_color_OkabeIto()

Some takeaways!

• Depending on the goals of our analysis, we may want to select a domain-specific methodology.
• We’ve reduced the dimensions of our data from 10s of 1000s of genes to 58 latent variables that we could test for differences between subtypes.

Other resources

• This website explains PCA visually.
  
• Wattenberg et al. (2016) discusses how to use t-SNE properly with great visuals. (The lessons apply to UMAP as well, with a broad substitution of the n_neighbors parameter for perplexity.)
• Nguyen & Holmes (2019) lay out guidelines on choosing dimensions reduction methods.
  
• Freytag (2019) is a nice explanation and comparison of many different dimensionality reduction techniques that you may encounter.
• Raj (2019) - another nice introduction to dimensionality reduction

Session Info

sessionInfo()

R version 4.0.2 (2020-06-22)
Platform: x86_64-pc-linux-gnu (64-bit)
Running under: Ubuntu 20.04 LTS

Matrix products: default
BLAS/LAPACK: /usr/lib/x86_64-linux-gnu/openblas-openmp/libopenblasp-r0.3.8.so

locale:
 [1] LC_CTYPE=en_US.UTF-8       LC_NUMERIC=C               LC_TIME=en_US.UTF-8        LC_COLLATE=en_US.UTF-8    
 [5] LC_MONETARY=en_US.UTF-8    LC_MESSAGES=C              LC_PAPER=en_US.UTF-8       LC_NAME=C                 
 [9] LC_ADDRESS=C               LC_TELEPHONE=C             LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C       

attached base packages:
[1] stats     graphics  grDevices utils     datasets  methods   base     

other attached packages:
 [1] Rtsne_0.15      uwot_0.1.8      Matrix_1.2-18   forcats_0.5.0   stringr_1.4.0   dplyr_1.0.0     purrr_0.3.4    
 [8] readr_1.3.1     tidyr_1.1.0     tibble_3.0.3    ggplot2_3.3.2   tidyverse_1.3.0

loaded via a namespace (and not attached):
 [1] tinytex_0.25       tidyselect_1.1.0   xfun_0.16          haven_2.3.1        lattice_0.20-41    colorspace_1.4-1  
 [7] vctrs_0.3.2        generics_0.0.2     htmltools_0.5.0    yaml_2.2.1         base64enc_0.1-3    blob_1.2.1        
[13] rlang_0.4.7        pillar_1.4.6       withr_2.2.0        glue_1.4.1         DBI_1.1.0          dbplyr_1.4.4      
[19] modelr_0.1.8       readxl_1.3.1       matrixStats_0.56.0 lifecycle_0.2.0    munsell_0.5.0      gtable_0.3.0      
[25] cellranger_1.1.0   rvest_0.3.6        evaluate_0.14      labeling_0.3       knitr_1.29         irlba_2.3.3       
[31] fansi_0.4.1        broom_0.7.0        Rcpp_1.0.5         scales_1.1.1       backports_1.1.8    jsonlite_1.7.0    
[37] RSpectra_0.16-0    FNN_1.1.3          farver_2.0.3       fs_1.5.0           hms_0.5.3          digest_0.6.25     
[43] stringi_1.4.6      grid_4.0.2         cli_2.0.2          colorblindr_0.1.0  tools_4.0.2        magrittr_1.5      
[49] crayon_1.3.4       pkgconfig_2.0.3    ellipsis_0.3.1     xml2_1.3.2         reprex_0.3.0       lubridate_1.7.9   
[55] assertthat_0.2.1   rmarkdown_2.3      httr_1.4.2         rstudioapi_0.11    R6_2.4.1           compiler_4.0.2